Perturb-Then-Diagonalize Vibrational Engine Exploiting Curvilinear Internal Coordinates

The present paper is devoted to the implementation and validation of a second-order perturbative approach to anharmonic vibrations, followed by variational treatment of strong couplings (GVPT2) based on curvilinear internal coordinates. The main difference with respect to the customary Cartesian-based formulation is that the kinetic energy operator is no longer diagonal, and has to be expanded as well, leading to additional terms which have to be taken into proper account. It is, however, possible to recast all the equations as well-defined generalizations of the corresponding Cartesian-based counterparts, thus achieving a remarkable simplification of the new implementation. Particular attention is paid to the treatment of Fermi resonances with significant number of test cases analyzed fully, validating the new implementation. The results obtained in this work confirm that curvilinear coordinates strongly reduce the strength of inter-mode couplings compared to their Cartesian counterparts. This increases the reliability of low-order perturbative treatments for semi-rigid molecules and paves the way toward the reliable representation of more flexible molecules where small- and large-amplitude motions can be safely decoupled and treated at different levels of theory.


S1 Kinetic energy operator in curvilinear coordinates
In the following, the mathematical steps to convert the standard equation of the internalbased kinetic energy operator T into a form where different contributions are separated, are shown.
For readability, the initial expression S1 is recalled here, where s = {1, . . . , s n } is a complete and non-redundant set of internal coordinates and G is the determinant of the Wilson G matrix.
In order to rewrite Eq. S1, the application of the differential operator T to a generic function Φ(s) must be explored: Now it is possible to notice that the fourth and fifth terms within parenthesis can be simplified, S-4 and upon rearrangement, the kinetic energy operator can be written as Finally, a further algebraic manipulations yields the following, reference expression, where V ′ is known as extra-potential term, whose expression is reported below:

S2 Wilson GF method
In analogy with the treatment in terms of Cartesian-based normal coordinates (CNCs), the harmonic approximation of vibrations can be employed, through the so-called Wilson GF method. S2 Within this method, the G matrix is approximated with its value at the reference geometry (usually the equilibrium one). As a result, the dependence of the G matrix on the internal coordinates is neglected, as well as that of its determinantG, so that, where the superscript "eq" indicates that the elements G ij are evaluated at the equilibrium geometry, and from here on will be omitted.
The treatment of potential energy is identical to that employed in terms of CNCs, the only difference being the set of coordinates used for the Taylor-series expansion, By assuming the condition of stationary point and shifting the origin of the PES to zero, the above equation truncated to the second order can be written as, where F is the Hessian matrix of the potential energy with respect to the internal coordinates. While the G matrix is easily obtainable on the basis of structural parameters (namely the molecular geometry and atomic masses), F can be calculated from the Cartesian Hessian matrix H x . In the most general case, also accounting for non-equilibrium geometries, where g s is the gradient of the potential expressed in internal coordinates, which can be calculated from its Cartesian counterpart g x , S3,S4 U is a (3N a × 3N a ) arbitrary matrix, with N a being the number of atoms. It is noteworthy that the contribution of translations and rotations can be factored out through the application of the projection matrix P = B † B. More specifically, g s and F can be obtained as, S4,S7 Once this procedure has been carried out, the harmonic vibrational Hamiltonian H (0) can be defined as, At variance with the expansion in terms of mass-weighted Cartesian displacements, the kinetic energy is not diagonal, the different coordinates being coupled by the G matrix. As a result, a set of coordinates diagonalizing both F and G, must be defined. For this purpose, a set of normal coordinates Q can be introduced, By inserting Eq. S15 in Eq. S14, and recasting the latter in a matrix form, we obtain the following expression, L is defined so that in this basis the F matrix becomes diagonal (Λ), while G is equal to the identity matrix. This corresponds to the resolution of the following equation, In other words, the calculation of the harmonic frequencies and internal-based normal coordinates (INCs) can be performed through the diagonalization of the GF matrix product.

S-7
It should be noted that for curvilinear coordinates the matrix to be diagonalized is not symmetric, implying that the normal coordinates do not form an orthogonal basis. However, as demonstrated by Myazawa, S8 Eq. S17 can be recast in a symmetric form characterized by the same eigenvalues, where the columns of the matrix G −1/2 L are the eigenvectors of the symmetric matrix G 1/2 FG 1/2 , and thus are orthogonal.
It is important to point out that the harmonic Hamiltonian H (0) introduced in Eq. S14, can be converted to the basis of the dimensionless normal coordinates q , leading to the following expression, where p i = −i∂/∂q i . Since the expression of H (0) through dimensionless units is equivalent to that obtained in the Cartesian-based formulation, the eigenvalues and eigenfunctions of H (0) are still given by Hermite polynomials multiplied by gaussian functions. Moreover, the second quantization formalism can be employed in the derivation of the equations stemming from the perturbative treatment.

S3 Conversion of energy units and Wilson G matrix derivatives
In this section, the expressions required to convert both energy and Wilson G matrix derivatives to wavenumbers is shortly discussed. For the purpose, it is necessary to introduce the so-called dimensionless normal coordinates q i (i = 1, . . . , N , with N being the number of normal modes) in terms of the mass-weighted normal coordinates Q i (in √ amu·Bohr), S-8 where γ i = 2πcω i /ℏ, c is the light speed (in cm/s), ω i is the harmonic wavenumber related to the i-th mode (in cm −1 ) and ℏ is the reduced Plank constant (in amu·Bohr 2 s −1 ). The conjugate moment of q i is p i = −i∂/∂q i .

S3.1 Kinetic-energy derivatives
Let us consider a generic element G ij of the Wilson G matrix and let us define its derivative with respect to the mass-weighted normal coordinates Let us now introduce the g matrix, obtained by converting the G matrix to wavenumbers, and its derivative with respect to the dimensionless normal coordinates: The transformation between the two derivatives is given by By considering the zero-th order term and the first-and second-order derivatives of Eq. S23, corresponding to the contributions appearing in the perturbative expansion of the kinetic energy operator up to the second-order, the following expression arises Furthermore, the Wilson G matrix is equal to the identity matrix at the equilibrium geometry:

S3.2 Potential-energy derivatives
Let us consider a generic anharmonic force constant with respect to the mass-weighted normal coordinates, The corresponding counterpart expressed in wavenumbers is The transformation between the two derivatives is given by By considering the second-, third-and fourth-order derivatives of Eq. S28, corresponding to the contributions appearing in the perturbative expansion of the potential energy up to the second-order, the following expressions are obtained: Let us stress that F ij = 4π 2 c 2 ω 2 i δ ij at the equilibrium geometry, so that

S4 χ matrix for Cartesian normal modes
The anharmonic χ matrix in the Cartesian-based formulation of VPT2 is The expressions reported above can be rearranged in order to isolate the potentially resonant contributions: S-11

S5 Resonant terms in the χ matrix
In this section, the potential (χ V ), kinetic (χ T ) and cross (χ × ) contributions to the internalbased formulation of the χ matrix, whose potentially resonant contributions have been isolated, are reported.
Such expressions are obtained starting from Eqs. 34 to 39 and applying the following identities,

S5.1 Potential term
The partial fraction decomposition of the potential term is given by,

S5.2 Kinetic term
The partial fraction decomposition of the kinetic term is given by, where s ijk = g ij,k + g ik,j + g jk,i r ijk = g ij,k − g ik,j − g jk,i

S5.3 Cross term
The partial fraction decomposition of the cross term is given by, g ii,k f jjk + f iik g jj,k ω k